function [duration, duration_nbr, divorcees_nbr, model_sahw_t]=moments_postestimation(n_g, T, lt, dt, H_max, disc, stoch_earn_f, stoch_earn_m, m_eta_f, m_eta_m, n_eta_f, n_eta_m, home_prod, match_qual, m_th, n_th, gamma, L_max, params, data_moments, data_f, data_m, W)

%----------------------- Preliminaries -----------------------------------%
parest(1,:)=params(1:9); 
parest(2,:)=params(10:18); 
sigma_f=[params(19) params(19) params(19) params(20) params(20) params(20) params(21) params(21) params(21)]; 
sigma_m=[params(22) params(23) params(24) params(22) params(23) params(24) params(22) params(23) params(24)]; 
match_qual(:,2)=params(25:33)'; 

%- Wage rate
a0m_1=params(34);
a0m_2=params(35);
a0m_3=params(36);

a0f_1=params(37);
a0f_2=params(38);
a0f_3=params(39);

%- Returns to experience
a1m_1=params(40);
a1m_2=params(41);
a1m_3=params(42);

a1f_1=params(43) ;
a1f_2=params(44) ;
a1f_3=params(45) ;

%- Returns to experience^2 and beyond
a2m_1=params(46);
a3m_1=0;

a2m_2=params(47);
a3m_2=0;

a2m_3=params(48);
a3m_3=0;

a2f_1=params(49);
a3f_1=0;

a2f_2=params(50);
a3f_2=0;

a2f_3=params(51);
a3f_3=0;

%- Returns to sahw
bm_1=params(52);
b2m_1=0;
bm_2=params(53);
b2m_2=0;
bm_3=params(54);

b2m_3=0;
bf_1=0;
b2f_1=0;
bf_2=0;
b2f_2=0;
bf_3=0;
b2f_3=0;

det_earn_m_1=[a0m_1, a1m_1, a2m_1, a3m_1, bm_1, b2m_1];
det_earn_m_2=[a0m_2, a1m_2, a2m_2, a3m_2, bm_2, b2m_2];
det_earn_m_3=[a0m_3, a1m_3, a2m_3, a3m_3, bm_3, b2m_3];

det_earn_m=[det_earn_m_1;det_earn_m_2;det_earn_m_3;det_earn_m_1;det_earn_m_2;det_earn_m_3;det_earn_m_1;det_earn_m_2;det_earn_m_3];

det_earn_f_1=[a0f_1, a1f_1, a2f_1, a3f_1, bf_1, b2f_1];
det_earn_f_2=[a0f_2, a1f_2, a2f_2, a3f_2, bf_2, b2f_2];
det_earn_f_3=[a0f_3, a1f_3, a2f_3, a3f_3, bf_3, b2f_3];

det_earn_f=[det_earn_f_1;det_earn_f_1;det_earn_f_1;det_earn_f_2;det_earn_f_2;det_earn_f_2;det_earn_f_3;det_earn_f_3;det_earn_f_3];

% Pareto Weights
pweight = [params(55:63); params(64:72); params(73:81); params(82:90)]; 

%------------------------------ Moments ----------------------------------%
model_behavioral_moments_bymkt = zeros(42, 4); 
U_mc_bymkt = zeros(n_g,4,4); 
U_mc_tot_bymkt = zeros(n_g,2,4); 

duration=zeros(n_g,4); 
duration_nbr=zeros(n_g,4); 
divorcees_nbr=zeros(n_g,4); 

model_sahw_t = zeros(n_g,4,4); 

I = 10000;

lcd_bymkt=zeros(n_g,I,T-lt+1,4); 

for mkt=1:4     
parest(3,:)=pweight(mkt,:); 

[dec, U_f, U_m, U_f0, U_m0, U_ftot, U_mtot] = solution_backwards(n_g, T, lt, dt, H_max, disc, det_earn_f, det_earn_m, stoch_earn_f, stoch_earn_m, m_eta_f, m_eta_m, n_eta_f, n_eta_m, parest(1,:), home_prod, parest(2,:), match_qual, m_th, n_th, gamma, L_max, parest(3,:), sigma_f, sigma_m);
  
hw = ones(n_g,I,T+1); 
lc = zeros(n_g,I,T+1); 
lcd=zeros(n_g,I,T-lt+1); 
hwd=zeros(n_g,I,T-lt+1); 

n_p=4; 

married_moments=zeros(n_g,n_p);
singles_moments=zeros(6,1); 

for g=1:n_g
married_moments(g,1)=(exp(U_f(g))/U_ftot(g))*data_f(g,mkt); 
married_moments(g,2)=(exp(U_m(g))/U_mtot(g))*data_m(g,mkt);
end

singles_moments(1,1)=(exp(U_f0(1))/U_ftot(1)); 
singles_moments(2,1)=(exp(U_f0(4))/U_ftot(4));
singles_moments(3,1)=(exp(U_f0(7))/U_ftot(7));

singles_moments(4,1)=(exp(U_m0(1))/U_mtot(1)); 
singles_moments(5,1)=(exp(U_m0(2))/U_mtot(2));
singles_moments(6,1)=(exp(U_m0(3))/U_mtot(3));

%- Set the seed

seed=1812;
rng(seed); 

wft_rnd_all=rand(3,I,T); 
wft_rnd_all_rep=[wft_rnd_all(1,:,:) ; wft_rnd_all(1,:,:) ;wft_rnd_all(1,:,:) ;wft_rnd_all(2,:,:) ; wft_rnd_all(2,:,:) ;wft_rnd_all(2,:,:); wft_rnd_all(3,:,:) ; wft_rnd_all(3,:,:) ;wft_rnd_all(3,:,:) ];

wmt_rnd_all=rand(3,I,T); 
wmt_rnd_all_rep=[wmt_rnd_all(1,:,:) ; wmt_rnd_all(2,:,:) ;wmt_rnd_all(3,:,:) ;wmt_rnd_all(1,:,:) ; wmt_rnd_all(2,:,:) ;wmt_rnd_all(3,:,:) ;wmt_rnd_all(1,:,:) ; wmt_rnd_all(2,:,:) ;wmt_rnd_all(3,:,:) ];

tht_rnd_all=rand(9,I,T); 

wft_rnd=zeros(I,T); 
wmt_rnd=zeros(I,T); 
tht_rnd=zeros(I,T); 

wft_gridpt=zeros(I,T);
wmt_gridpt=zeros(I,T);
th_gridpt=zeros(I,T);

w_node=zeros(n_g,n_eta_f);
m_node=zeros(n_g,n_eta_m);

for g=1:n_g 

%-Female wages
[wf_node(g,:),~,cdf_fu,cdf_fc]=wageproc_fn(stoch_earn_f(g,:), m_eta_f, n_eta_f);

%-Male wages
[wm_node(g,:),~,cdf_mu,cdf_mc]=wageproc_fn(stoch_earn_m(g,:), m_eta_m, n_eta_m);

%-Theta process
[~,~,cdf_thu,cdf_thc]=thetaproc(match_qual(g,:), m_th, n_th);

wft_rnd(:,:)=wft_rnd_all_rep(g,:,:); 
wmt_rnd(:,:)=wmt_rnd_all_rep(g,:,:); 
tht_rnd(:,:)=tht_rnd_all(g,:,:);
%- Female wage nodes

for i=1:I   
       for t=lt:T 
 
           if t==lt; 
               
        for f=1:n_eta_f 
            if wft_rnd(i,t)<=(cdf_fu(1,f));
               wft_gridpt(i,t)=f;
               pastnode=f;
               break
            end
        end
  
           else 
               
        for f=1:n_eta_f
            if wft_rnd(i,t)<= (cdf_fc(pastnode,f));
            wft_gridpt(i,t)=f;
            newnode=f;
            break
            end
        end
        
 
        pastnode=newnode;
                   
           end
           
       end
end

%- Male wage nodes

for i=1:I   
       for t=lt:T 
             
           if t==lt;
               
        for m=1:n_eta_m 
            if wmt_rnd(i,t)<=(cdf_mu(1,m));
               wmt_gridpt(i,t)=m;
               pastnode=m;
               break
            end
        end
        
            else
               
        for m=1:n_eta_m 
            if wmt_rnd(i,t)<= (cdf_mc(pastnode,m));
            wmt_gridpt(i,t)=m;
            newnode=m;
            break
            end
        end
  
        pastnode=newnode; 
        
           end
           
       end
end


%- Theta nodes
for i=1:I   
       for t=lt:T 
           
           if t==lt;
        for th=1:n_th 
            if tht_rnd(i,t)<=(cdf_thu(1,th));
               th_gridpt(i,t)=th;
               pastnode=th;
               break
            end
        end
 
           else
               
        for th=1:n_th  
            if tht_rnd(i,t)<= (cdf_thc(pastnode,th));
            th_gridpt(i,t)=th;
            newnode=th;
            break
            end
        end

        pastnode=newnode; 
        
           end
           
       end
end

hw(g,:,lt+1:T+1) = 0; 
for i=1:I
   
    for t=lt 
        lc(g,i,t) = dec(g,hw(g,i,t), wft_gridpt(i,t), wmt_gridpt(i,t), th_gridpt(i,t), t);
        if lc(g,i,t) == 2
            hw(g,i,t+1) = hw(g,i,t) + 1; 
        elseif lc(g,i,t) == 1
            hw(g,i,t+1) = hw(g,i,t);
        end
    end
    
    for t=lt+1:T 
        
        if lc(g,i,t-1)==3 
            lc(g,i,t)=-1;            
        elseif lc(g,i,t-1)==-1 
            lc(g,i,t)=-1;
            
        else 
        lc(g,i,t) = dec(g, hw(g,i,t), wft_gridpt(i,t), wmt_gridpt(i,t), th_gridpt(i,t), t); 
        end                                                                   
        if lc(g,i,t) == 2
            hw(g,i,t+1) = hw(g,i,t) + 1;
        elseif lc(g,i,t) == 1
            hw(g,i,t+1) = hw(g,i,t);
        else
            hw(g,i,t+1) = 0; 
        end
        
    end
end

lcd(g,:,:)=lc(g,:,lt:T); 
hwd(g,:,:)=hw(g,:,lt:T);

lcd_bymkt(g,:,:,mkt)=lcd(g,:,:);

wft_gridpt_bymkt(g,:,:,mkt)=wft_gridpt(:,:); 
wmt_gridpt_bymkt(g,:,:,mkt)=wmt_gridpt(:,:);

fracd=lcd==3;
frachw=lcd==2;
divorced=(lcd==3 | lcd==-1);

ever_divorced(:,:)=divorced(g,:,:);
ever_divorced_indx=sum(ever_divorced,2);
still_married(:,:)=(lcd(g,:,:)==1 | lcd(g,:,:)==2); 
still_married(ever_divorced_indx==0,:)=[]; 
nbr_div=size(still_married,1);

married_moments(g,3)=sum(sum(frachw(g,:,1:T)))/(I*T-sum(sum(divorced(g,:,1:T)))); 
married_moments(g,4)=sum(sum(fracd(g,:,1:T)))/I; 

duration(g,mkt)=sum(sum(still_married(:,1:T)))/nbr_div; 
duration_nbr(g,mkt)=sum(sum(still_married(:,1:T)));
divorcees_nbr(g,mkt)=nbr_div; 

clear still_married

model_sahw_t(g,1,mkt)=sum(sum(frachw(g,:,1:2)))/(I*2-sum(sum(divorced(g,:,1:2))));   
model_sahw_t(g,2,mkt)=sum(sum(frachw(g,:,3:4)))/(I*2-sum(sum(divorced(g,:,3:4))));   
model_sahw_t(g,3,mkt)=sum(sum(frachw(g,:,5:6)))/(I*2-sum(sum(divorced(g,:,5:6))));  
model_sahw_t(g,4,mkt)=sum(sum(frachw(g,:,7:10)))/(I*4-sum(sum(divorced(g,:,7:10))));  

end 

end 

end




